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A Method for the Numerical Integration of Differ- 
ential Equations of Second Order Without 
Explicit First Derivatives * 

Rene de Vogelaere 2 



This is a fourth-order step-by-step method based on difference formulas. The case of 
a single equation is discussed in detail. The use of this method for automatic computer is 
considered. 

1. Introduction 

This is a fourth-order step-by-step method based on difference formulas. The start of 
the computation is easy. Only the case of a single differential equation 

g=/(*,?y) (i) 

will be considered. The generalization to a set of differential equations 

d 2 y 



dx 



j=fi(x,Vk) i,k=l,2, . . . ,n 
is immediate. Application to electronic computing is discussed. 

2. General Formulas 

We first note some general formulas that may be deduced, for instance, from the expres- 
sions (1.32), (1.38), (1.39), and (4.15) of L. Collatz [l]: 3 



f f V (x)dx 2 = ™j(x) - (2) /(0) - <»/(0)s 
Jo Jo 



4 ! Eft,,WA'i- f (2) 

=h 2 J2(-\yP 2 , p (-u)A*f (3) 

P=0 

where x=uh, 

P 2 ,o(u)=u 2 /2; P 2 ,i(u)=u*/Q) P 2 . 2 (u)=(2<a 3 + u')/24; P 2 .3N = (20u 3 +15u 4 +3i/ 5 )/360; ... (4) 

and 

(n)f = T (n-i )f(lx (o)y=y ? n=l } 2. 



Finally, 



(1)/ CD/. — Af9/._J-I A2/ L A4 



>/ 2 -a>/ =/ l (2/ 1 + ±A 2 /„-^/_ 1 + . . .)• (5) 



i Work performed in part when the author was a guest worker at the National Bureau of Standards, November 1953, and in part at the Uni- 
versity of Louvain, Belgium, 194S, toward the author's Dissertation. 

2 Department of Mathematics, University of Notre Dame, lnd. 

3 Figures in brackets indicate the literature references at the end of this paper. 
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The degree of approximation of the formulas given below is written in terms of the interval 
h; the symbol [n] is equivalent to the usual o(h n ). The following remarks are needed: 

(1) /, computed from its definition, is of the same order of approximation as that of the y 
used ; 

(2) if A p+ \f is neglected in (2) or (3), these expressions are correct to the (p + 2)nd order; 

(3) if A 4 /_! is neglected in (5), this formula is correct to the fourth order. In the notation 
introduced, the differential equation under consideration can be written as 

we shall also write 

dx J 

3. Particular Formulas 

When x=h, and if A 2 /_ 2 is neglected, (2) becomes 

y 1 =yo+hz +h 2 (-f_ l +4f )/^, {3} (6) 

provided y { 3}, z G {2}, and/ {l} ? /-i{l}- 

If A/_i is neglected in the same formula (2), we have 

yi=yo+hzo+h 2 f /2. {2} (7) 

When x=h, and if A 2 / is neglected, (3) becomes 

yi=yo+hzo+h\2f +fi)/^ {3} (8) 

provided y { 3}, z {2], and/ {l},/i{l}- 

The same formula yields, for x=2h, and with A 3 / neglected, since P 2 ,2 (—2) =0, 

y 2 =y + 2hz +h 2 (2f +±f 1 + 0f 2 )/3, {4} (9) 

and provided ?/o{4}, z {3} , and/ {2},/i{2}. 

Neglecting A 4 /_i, (5) gives the Simpson formula 

^=^o+AC/o+4/ 1 +/ 2 )/3 {4} (10) 

if o {4}and/o{3},/i{3},/ 2 {3}.' < 

Finally, when x=2ht, and if A 3 / is neglected, (3) becomes 

y=y +2hz Q t+2h%f+W(-3fo+\fi-f2)t*+W {4} (11) 

4. Integration Method 

We are now ready to explain our method for solving differential equations of second order 
that do not contain the first derivatives explicitly. The starting value of the independent 
variable will be taken as zero. Each step is subdivided into two intervals of length h. 

The general step will be considered first. Therefore, we will suppose that /=/_ a { 3} is 
known for x=(2n— l)k, and that y=y {4},f=f {4:}, z=z {4:}, are known for 2=2nh. 

Then for x=(2n+l)h, y=yi{3} is given by (6) and/=/i{3} is given by (1), whereas for 
x=(2n+2)h, y=y 2 {±} is given by (9),/=/ 2 {4} by (1) and z=z 2 {4} by (10). 

These values will then serve for the next step. The process is of fourth order for the 
function as well as for its first derivative. 

Two equivalent methods may be used to compute the first step. For x=0, and with 
initial conditions y = y and z=z Q , f=fo is computed by means of the differential equation. 
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Then 

(a). A preliminary calculation may be done with (7) and (1) so as to give 

y=y 1 {2} and/=/ 1 {2},fora?=A. 

Then formula (8) may be used with/i=/ij this gives 

y=yi{3} and/=/i{3} ; for x=h. 

The computation of the values for x=2h is the same as for the general step (2n J r2)h. 

(b). Or, alternatively, it is also possible first to go backwards with the formula equivalent 
to (7), 

y-i=y<>-hz +h 2 f /2. {2} 

The process for the general step may then be used for x=h and x=2h because (6) requires 
only a knowledge of /_ t to the first order and this function is known to the second. 

Table 1 gives the coefficients of the z and/ functions, (the coefficient of y is always one) 
for certain intervals between 0.50 and 0.06; if the interval to be chosen is outside this range, 
it is preferable, especially with the use of desk machines, to change the independent variable; 
this will prevent rounding off errors or the need of computation with an excess of decimal 
places. 

Table 1. 
For the users' convenience, certain entries have been repeated. 















Coefficients in formulas— 














Inter- 
vals 


(6) 


(S>) 


(10) 


(7) 


(8) 


(7) 




zo 


-/-. 


H 


20 


h 


/, 


/o+4/i+/2 


20 


h 


20 


/o 


/. 


-20 


/o 


h 


h 


fc2/6 


2/J-73 


2ft 


2/12/3 


4/^2/3 


fc/3 


h 


fc2/2 


h 


/i.2/3 


A 2/6 


h 


A2/2 


0. 060000 


0. 060000 


0. 000600 


0. 002400 


0. 120000 


0. 002400 


0. 004800 


0. 020000 


0. 060000 


0.001800 


0. 060000 


0.0012(H) 


0. 000600 


0. 060000 


0. 001800 


. 075000 


. 076000 


. 000938 


. 003750 


. 150000 


. 003750 


. 007500 


. 025000 


. 075000 


.002813 


.075000 


.001875 


.000938 


. 075000 


.002813 


. 090000 


. 090000 


.001350 


. 005400 


. 180000 


. 005400 


. 010800 


.030000 


.090000 


.()( 14050 


.090000 


. 002700 


. 001360 


. 090000 


.004060 


. 100000 


. 100000 


.001007 


. 006667 


. 200000 


. 006667 


. 013333 


. 033333 


. 100000 


.005000 


. 100000 


. 003333 


.001667 


. 100000 


. 005000 


. 120000 


.120000 


. 002400 


. 009600 


. 240000 


. 009600 


.019200 


. 040000 


. 120000 


. 007200 


. 120000 


. 004800 


. 002400 


. 120000 


. 007200 


. 125000 


. 125000 


. 002604 


. 010417 


. 250000 


.010417 


. 020833 


.041667i 


.125000 


. 007813 


. 125000 


. 005208 


002001 


. 125000 


. 007813 


. 150000 


. 150000 


. 003750 


. 015000 


. 300000 


. 015000 


. 030000 


. 050000 


. 1.50000 


.011250 


. 150000 


. 007500 


. 003750 


. 150000 


.011250 


. 180000 


. 180000 


. 005400 


. 021600 


. 360000 


. 021600 


. 043200 


. 060000 


. 180000 


.016200 


. 180000 


. 010800 


. 005400 


. 180000 


. 016200 


. 200000 


. 200000 


. 006667 


. 026667 


. 400000 


. 026667 


. 053333 


. 066667 


. 200000 


. 020000 


. 200000 


.013333 


. 006667 


. 200000 


. 020000 


. 240000 


. 240000 


. 009600 


. 038400 


. 480000 


. 038400 


. 076800 


. 080000 


. 240000 


. 028800 


. 240000 


.019200 


. 009600 


. 240000 


. 028800 


. 250000 


. 250000 


. 010417 


. 041667 


. 500000 


. 041667 


. 083333 


. 083333 


. 250000 


. 031250 


. 250000 


. 020833 


.010417 


. 250000 


. 031250 


. 300000 


. 300000 


. 015000 


. 060000 


600000 


. 060000 


. 120000 


.100000 


. 300000 


. 045000 


. 300000 


. 030000 


. 015000 


. 300000 


. 045000 


. 360000 


. 360000 


. 021600 


. 086400 


. 720000 


. 086400 


. 172800 


. 120000 


. 360000 


. 064800 


. 360000 


. 043200 


. 021600 


. 360000 


. 004800 


. 400000 


. 400000 


. 026667 


. 106667 


. 800000 


. 106667 


. 213333 


. 133333 


. 400000 


. 080000 


. 400000 


. 053333 


. 026667 


. 400000 


.080000 


. 420000 


. 420000 


. 029400 


. 117600 


. 840000 


. 117600 


. 235200 


. 140000 


. 420000 


. 088200 


. 420000 


. 058800 


. 029400 


. 420000 


. 088200 


. 480000 


. 480000 


. 038400 


. 153600 


. 900000 


. 153600 


. 307200 


. 160000 


. 480000 


. 115200 


. 480000 


. 076800 


. 038400 


. 480000 


. 115200 


. 500000 


. 500000 


. 041667 


. 166667 


1.000000 


. 166667 


. 333333 


. 1666671 


. 500000 


. 125000 


. 500000 


. 083333 


.041667 


. 500000 


. 125000 



5. Application to High-Speed Computors 
For this application, it is more convenient to introduce the quantities 



hz 2 



J 2ni 



-^h 2 f 2n —F2n) 



-Mhn 



+ 1- 



2»+l- 



From our formulas it is easy to deduce that only four cells are needed to store the values of y, 
Z, and F. We shall give here a possible arrangement of the calculation so as to have only one 
multiplication when the jF's are computed. 
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Thus for the general step, if in the memory cells, P, Q, R, S, we have at a certain moment 
the following situation 4 

CP)=2/o, (Q)=Zo, (B)=F , (S)=F_ U 
the following five orders 

(Q) + (B)^Q, (P) + (<2)-P, \(S)-*S, (R)-(S)^S, (P) + (S)^S, 

put ?/i in S or in the y cell of the subroutine corresponding to (1). This subroutine produces- 
Fij which is stored in S. The two orders 

(Q) + (S)-*Q, (P) + (Q)^P, 

then put y 2 in P, and the subroutine corresponding to (1) computes F 2 , which is then stored in 
R, Finally, the order 

(Q) + (P)^Q, 

puts Z 2 in Q. 

To start the computation the method (b) is preferable because once the backward step 

3 

y_i=y — Z + -F 

has been made, we immediately enter the general routine. 

6. Change of Interval 

Because /_i needs only to be of first order in (6) , there is no difficulty in changing the 
interval from hi to h 2 for the value x of the independent variable. For, if j_ x denotes the value 
of /for the last step (x—x — hi) and/_i, the value of /which is needed (#— £ — A 2 ), we have 

/-i""/o - J-i"/o r 1 \ 

lh h 2 {l] 

Inserting this in (6) we have 

y^yo+h^-f (^)JZH-| W (l+l g/o. {3} (12> 

In this form, the quantities in the parenthesis are multipliers of the coefficients of (6) for a new 
interval h 2 . 

The formula (8) may be deduced from (12) with the particular value hi=—h 2 and with 
/_i replaced by/i. 

7. Error of the Process 

Only the error of one general step will be considered. Moreover, since the first derivative 
has been calculated to the same order as the function itself, i. e., the fourth, only the fifth-order 
term of the function is of interest. If y , Zo, and/ are correct, this term is given by the first 
term neglected when (9) is derived, namely, 

^ 2 A 3 /o, (13) 

since /1 is known to the third order in h, regardless of which of the formulas, (6) or (8) is used 



* (X) means contents of cell X. 
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to compute y x . An estimate of the error is thus given by 



2_ 

45 



KSHKS)/ 



However, it is (13) that provides a practical method of estimating the total error. 
To check the computation, two methods may be used. 

(a) When a first integral of the equation or of the set of equations is known, this may be 
used at each step as a check because z has been computed. But confidence in this as a check 
should not be too great. Although, if the integration is part of a larger program, this check will 
be in general sufficient. 

(b) If an isolated integration has to be done, greater care may be needed. The following 
procedure is then proposed. Recompute y x (which is known with third-order accuracy) by 
means of the interpolation formula (11) considered backward, so as to use y 2 , 2 2 , and/ 2 - 

y *=y 2 - hz2+ ± /,2 ( _y o+6/i+7 y 2)> (15) 

The difference is, for the general step, approximately the first neglected term in (6), 

and should be small and vary slowly. 

The accuracy of the calculation may also be deduced from the comparison of the differ- 
ences of yf—yi for two consecutive steps, 

(yf-y 1 h n+ 2-(yf-yi)2 n =lh i .2h(^^ n>o, (17) 

which is 8/45, the allowed error. The variation of yf—yi from one step to the other must 
therefore remain smaller than the 45/8 of the allowed error for each step. 

8. Comparison With Other Methods 

It is not our intention to make a detailed comparison here with similar methods — ordinary 
step by step method, Gauss-Jackson method, Runge-Kutta method. This task is very delicate, 
and it is necessary to take care of the higher order terms neglected, the precision needed, and 
the purpose of the integration. 

But we should like to draw attention to the fact, that in most of the other step-by-step 
methods, predictor and corrector formulas are used for the same value of the independent 
variable; here, a formula of the corrector type (9) permits us to advance. Moreover, the 
advantages are evident — it is easy to start the integration, and there is no problem in changing 
the interval. 

However, since it is possible to make a comparison with the Runge-Kutta method, we 
should like to develop this point a little more. It will be seen that the proposed method pos- 
sesses some of the advantages of the Runge-Kutta method without having some of its incon- 
veniences (see, for instance, Milne [4]). The Runge-Kutta method may be characterized by 
the fact that each step is independent and that functions and derivatives are calculated in 
general four times for the fourth-order process. Also, in contrast to the step-by-step processes, 
a maximum accuracy is not sought for each computation, the process being based on compensa- 
tion, at the end, of all errors made during the step. But for the case of differential equations of 
second order of the type at present under consideration, there is one set of formulas for which 
only three calculations are needed (see Collatz, p. 33-34), and for which, by accident, the 
function but not the first derivative is calculated each time with maximum accuracy. In this 
case the error obtained by comparison with the value y 2 of the Taylor expansion may be written 
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in the form 

yt-y*= A A 5 |Vf 6 (|£)/>] /, (i 8) 

as may be deduced from the two formulas following (3.13) of reference [1]. In these formulas 
we set with n=2, v=0, q=2, and replace h by 2h. 

In the case of a single differential equation of type (1) the operators jPand D are defined by 

HI + -£)X).- 

For the proposed method, (14) may be written with the same notation 

Having made this remark, we may now point out that the new method also subdivides the 
process in steps. These steps are not independent, but the dependence is small (first order 
only on the second derivative), and because of this small dependence, only two computations 
are needed instead of three without making the change of interval difficult. 

It appears also from a comparison of (18) and (19) that the function is better approximated 
in general, than with the Runge-Kutta method. This is not the case for the first derivative, 
which is determined in the Runge-Kutta method with fifth-order accuracy; however, this fifth- 
order accuracy may be completely illusory. 

Finally, we should like to point out that the fourth-order interpolation formula (11) bears 
a great resemblance, to an analogous formula developed by Lemaitre for the Runge-Kutta 
methods [3]. 

9. Example 

This method has been used principally for differential equations arising in dynamical prob- 
lems with two degrees of freedom. Therefore, we shall give an example of this kind that is 
taken from the problem of primary cosmic rays (cf . references of [2]) . 

The equations written with the above notations are 
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?2 

-j2=(e~ 2vi cos2 2/2— 1— tan 2 y 2 ) tan y 2 , 



and the first integral is 



( c ^) 2+ ( d dx) 2=ae2vi ~ 1 ~ tan2 y*+ 2e ~ Vl - e ~ 2Vi cos2 2/2. 

Table 2 presents the computation for the initial conditions 

y lt0 = .448080, 2/2.0=0.000000, y [ =0.000000, ?/2,o = 0.206279. 

The value a - = 0.070598 is deduced from the first integral. 

The method (a) of section 4 was used for starting. The results in A were obtained by 
using A=0.2 and retaining six decimal places. The results in B were obtained by using 
h=0A and retaining five decimal places. 
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Table 2. 



Case A 



X 


Vi 


V2 


/. 


h 


Z\ 


22 


Cx 


C 2 


Ei 


Ei 


S 


0.0 


0. 448080o 
. 446925 
. 446919 

. 443413 5 
. 437469 

. 428869 7 
. 417462 

. 4029632 
. 385128 

. 363697e 


0. OOOOOOo 
.041256 
. 041093 

. 081210b 
. 119383 

. L54666g 
. 186141 

. 2129969 
. 234476 

. 250095: 


-0.057742 
-.058631 
-. 058627 

-.061221 

-.065341) 

-.070734 

-. 076969 

-.083571 

-.089941 

-.095417 


-0.000000 
-.024492 
-.024394 

-.048620 
-.072382 

-.095146 
-. 116085 

-.134105 
-. 147896 

-. 156251 


0. OOJOOOo 


0. 206279i 













.2 












.2 






3 











.4 


-. 023564 7 


. 196532 7 











.6 


-14 


2 




.8 


-. 0497882 


. 1676463 


1 


1 





1.0 


-6 


8 




1.2 


-. 0806004 


. 1214067 


1 


2 





1.4 





19 




1.6 


-.1165174 


. 0626104 


1 


2 












Case B 



C, 



Ci 



Fa 



Ez 



OA) 
.4 
.4 



1.6 

2.0 



2.4 

2.8 



0. 44808:* 
14346 
1 1337 

12885 9 
. 40303 

. 36366s 
. 30912 

. 23888« 
. 15379 

. 06010b 



O.OOOI 

. 08251 

.MSll'.t 

. 15465, 
. 21287 

. 25005 2 
. 26204 

. 25096; 

.2 11)37 

. L7624o 



-0.05771 
-.0(1130 
-.0(1123 

-.07073 
-.08354 

-. 09541 
-. 10036 

-.09175 
-.06102 



-0.00000 

-.04943 

.04861 

-.09514 
-. 13402 

. L5621 
-.15266 



O.OOOOOn 



04978a 

. 11(11% 



. l(17(17i 
.062682 



-. 12404 

-. 07788 



. 19496g 
-.03068 -.24009] 



.056103 

.11826c 



2 

(i 



-12 



4 
13 



C\ and C2 are the fourth-order verification terms y£i— 2/1,1 and y$ t \— 2/2,1, computed by (16). 
E { and E 2 are the estimates of error deduced from (13). S is the error in the first integral. 
The quantities C h C2, E u E 2j and £ are to be multiplied by 10~ 6 in case A and by 10~ 5 in 
caseB. 

10. Remarks 

(1). It should be noted that (17) is not true for n=0, i. e., lor the first two strps; for, 
instead of (6), (8) was used, for which 

By comparison with (16) it appears that the first C is about minus one-third times the second 
one. 

(2). One extra figure is used fovy it 2n and 2 it 2 n (i=l,2) to lessen the rounding-off error: 
this changes very little in the computation, except that the second term in (9) must be written 



/2h 



) (10*o). 



\10> 
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